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Abstract 

A system of autonomous differential equations with a stable limit cycle and perturbed by small white 
noise is analyzed in this work. In the vicinity of the limit cycle of the unperturbed deterministic system, 
we define, construct, and analyze the Poincare map of the randomly perturbed periodic motion. We show 
that the time of the first exit from a small neighborhood of the fixed point of the map, which corresponds 
to the unperturbed periodic orbit, is well approximated by the geometric distribution. The parameter 
of the geometric distribution tends zero together with the noise intensity. Therefore, our result can be 
interpreted as an estimate of stability of periodic motion to random perturbations. 

In addition, we show that the geometric distribution of the first exit times translates into statistical 
properties of solutions of important differential equation models in applications. To this end, we demon- 
strate three examples from mathematical neuroscience featuring complex oscillatory patterns charac- 
terized by the geometric distribution. We show that in each of these models the statistical properties 
of emerging oscillations are fully explained by the general properties of randomly perturbed periodic 
motions identified in this paper. 

Keywords: Poincare map, random perturbations, limit cycle. 



1 Introduction 



Accurate description of many important dynamical phenomena in science and engineering is impossible 
without taking into account random factors. Even small random perturbations can transform deterministic 
dynamics in unexpected ways and create new asymptotic regimes, which are not present in the unperturbed 
deterministic system. Examples include large deviation type mechanisms of regular dynamics in randomly 
perturbed systems ifTTl l25l l39ll . stochastic resonance Q HI, stochastic stabilization ll24l l33l . and noise- 
induced synchronization ED . to name a few. Mathematical analysis of these and other related phenomena 
requires effective geometric theory of randomly perturbed dynamical systems, which belongs to the interface 
between two mathematical disciplines: the theories of dynamical systems and stochastic processes. The 
goal of this work is to extend the Poincare map method, the main geometric tool for studying stability of a 
periodic motion in deterministic systems to solutions of randomly perturbed differential equations. 
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The mathematical analysis of effects of random perturbations on dynamics of nonlinear systems was 
initiated by Pontryagin, Andronov, and Vitt in their pioneering paper PTTl . The first systematic investigation 
of stability of solutions of stochastic differential equations was undertaken by Khasminsky, who extended 
many methods of classical theory of ordinary differential equations (cf. PTl l23l ) to randomly perturbed 
systems ll24ll . Freidlin and Wentzell developed the asymptotic method of analysis of randomly perturbed 
dynamical systems based on large deviations estimates |[T8l . Asymptotic properties of solutions of stochastic 
ordinary differential equations were studied in ||43l [T9l . More recent approaches for studying randomly 
perturbed dynamical systems are based on ideas from the theory of dissipative dynamical systems [4] and 
those from the geometric theory for slow-fast systems [8]. A survey of asymptotic methods for randomly 
perturbed dynamical systems with a variety of applications is available in PBl . 

In qualitative theory of nonlinear dynamical systems, local stability analysis of invariant sets such as 
equilibria and periodic orbits plays an important role. There are many effective analytical techniques for 
studying stability of solutions of deterministic differential equations ll23l . For randomly perturbed systems, 
stability of invariant sets is reflected in the statistics of the times of the first exit from the corresponding 
domains. For domains containing a stable fixed point, the Freidlin-Wentzell theory of large deviations 
characterizes the asymptotics of the first exit time and the geometric location of the point of exit a random 
trajectory from the domain [18]. Furthermore, it is known that the limiting distribution of the first exit time 
is exponential |fl~3Tl . 

In the hierarchy of invariant sets of autonomous differential equations, equilibria are followed by pe- 
riodic orbits. The main tool for studying stability of a periodic orbit is the Poincare map (PM). The PM 
captures the behavior of trajectories in a typically small neighborhood of the periodic orbit. The fixed point 
of the Poincare map corresponds to the periodic orbit. Stability of this fixed point of the map translates 
into the stability of the periodic orbit. In this work, we consider an autonomous system of differential equa- 
tions in W 1 possessing a limit cycle and perturbed by state-dependent white noise process. Under general 
assumptions on the limit cycle, we derive the PM for the randomly perturbed system and study its stability. 
If the periodic orbit of the deterministic system is asymptotically stable, we show that the distribution of 
the first exit times is approximately geometric with the parameter of the geometric distribution tending to 
zero together with the noise intensity. This result can be interpreted as a form of stability of motion near the 
periodic orbit: while the trajectories of the PM eventually leave a small neighborhood of the fixed point, for 
small noise they remain in this neighborhood for a very long time with large probability. 

The geometric distribution of the first exit times resulting from stochastic perturbations of a stable limit 
cycle is important for understanding statistical properties of nonlinear oscillations generated by the randomly 
perturbed models. In fact, this work was motivated by our earlier analysis of a class of neuronal models in 
l25l . In conclusion of this paper, we discuss the applications of our results to irregular bursting and mixed- 
mode oscillations, two important oscillatory regimes encountered in conductance based models of neurons. 

The organization of the paper is as follows. Section [2] contains the formulation of the problem, pre- 
liminaries about the PM, and the statement of the main result. In Section [3] in a small neighborhood of the 
periodic orbit of the deterministic system, we derive a PM for the randomly perturbed problem. The analysis 
follows the derivation of the PM for deterministic ordinary differential equations ll23ll using the asymptotic 
expansions of solutions of stochastic differential equations on finite time intervals |9l [lOj QjO. Section [4] 
presents the analysis of the linearization of the PM. Here, we use the results of Kesten on the iterations of 
linear functions of random matrices (cf. [28 ]) to show that the first exit time of the linearized problem has 
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asymptotically geometric distribution. The analysis of this section is a substantial generalization of our pre- 
vious work f25t . In Section|5J we estimate the contribution of the nonlinear terms of the PM on the statistics 



of the first exit times. This concludes the proof our main result - Theorem |2.7| In Section|6j we illustrate the 
analysis of the previous sections with applications to three problems in mathematical neuroscience. To this 
end, we use two conductance -based models of single neurons and a model of electrically coupled network 
of pancreatic j3— cells. In the absence of noise, all these models exhibit stable periodic oscillations. Adding 
noise to these models results in more complex stochastic oscillatory regimes: irregular bursting and mixed- 
mode oscillations. We show that despite different mathematical formulations of these models and different 
forms of resultant oscillations, emergent stochastic oscillatory patterns in these models are formed due to the 
random perturbations of stable limit cycle oscillations. We verify numerically that the number of spikes in 
one burst and the number of small oscillations between consecutive spikes in the time series corresponding 



to these models are distributed approximately geometrically in accord with Theorem 2.7 We conclude this 



paper with a brief discussion of results of this work as well as related results in the literature in Section[7] 

2 Assumptions and results 
2.1 The model 

Consider an ordinary differential equation for x : M — > 

x = f(x), (2.1) 

where / : M. d+1 — > is twice continuously differentiable function, and all its partial derivatives up to 
the second order are uniformly bounded in Suppose 

x = u(t), u(t + 1) = u(t), t G R (2.2) 



is a nonconstant periodic solution of (2.1 1 with the least period 1 and nonvanishing derivative 

u(t) ^ 0, t G R. (2.3) 

The corresponding orbit is denoted by 

= {x = u(6), 6 G S 1 := R x /Z}. (2.4) 



Along with (2.1 ) we consider a randomly perturbed system 

±t = f(xt) + aP(x t )W t , (2.5) 

where Wt is a standard Wiener process in M^" 1 " 1 . Matrix P(x) G R( a!+1 ) >< ( rf+1 ) is nondegenerate for any 
x G The entries of P are continuously differentiable functions and all their partial derivatives are 



uniformly bounded in The noise intensity a > is considered a small parameter. Equation (2.5 I is 

understood in the sense of Ito BOl . 
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2.2 The local coordinates 



By Theorem VI. 1.1 of [23], in a small neighborhood of O, there is an orthonormal moving coordinate frame 



{v(9),z 1 (9),z 2 (9),...,z n _ 1 (9)}, v(9) 



\m\ 



9 £ S 1 . 



(2.6) 



such that the first vector v{9) points in the tangential direction to O. Denote 

Z{9) = co1(*l(0), . . . , z n _i(0)) € R^ 1 )** 

then 

x = u{0) + Z(9)p, (2.7) 

defines a smooth invertible transformation x ^ (9, p) £ 8 1 x R d in a small neighborhood of (cf . Theorem 
VI. 1.1 of (23). 



We use the moving coordinates to rewrite (2. 1 ) in a form more amenable for analysis. 



Lemma 2.1. In new coordinates (2.7), near O (2.5) has the following form 

9 t = l + a(9t) T p + 0(\ Pt \ 2 ) + a(h(9 t ) T + 0(\ Pt \))W ; . 
Pt = 



where 



a T (9) 

R(9) 
h(9) 

H(9) 



R(9 t )pt + 0{\p t \ z ) + a ( H{9 t ) 1 + 0(|p t |) ) W t , 



(Df(u(9))) s Z(9), 

i(9))Z(9) - Z T (9). 
- 2 P (u{9)Y f (u(9)) e 



\f(u(9))\ 

Z T (9)Df(u(9))Z(9) - Z T (9)Z'(9), 
1 



\f«0)) 
P(u(9)) T Z{9) € R( d+1 ) xd . 



(2.8) 
(2.9) 

(2.10) 

(2.11) 
(2.12) 

(2.13) 



where M s := 2 1 (M + M T ) stands for symmetric part of matrix M. 

Remark 2.2. By setting a = on the right hand side of (2.10) and ( |2.1 1| ), we obtain the deterministic 



equation (2. 1 1 rewritten in the moving coordinate frame near O. 



2.3 The Poincare map 

Our next goal is to construct the PM near the periodic orbit O. We first review properties of the PM of the 



deterministic system (2.1 1 and then turn to the PM for the randomly perturbed system (2.5 1. 
Throughout this subsection, we will use 



S ~ 5 = {x = «(0) + Z(0)p : M<*}, 



(2.14) 
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to denote a local section to periodic orbit O. Here, 5 > is assumed to be sufficiently small. 



Using the continuity of solutions of (2.1 ) with respect to initial data, for sufficiently small Si > 0, one 
can find 6% > Si such that for any initial condition xq := (0, po) £ Sq^, the first return time 



T(x ) = inf{t > : x(t) G 5 ,« a } 
is finite. The PM P : Sofa — > Sq : $ 2 is defined as follows 

P(/»o) : = Pr, 

where p {0jT} = Z(0) -1 X{ ,T}- 

The following properties of P are well-known (see, e.g., |42|): 



(2.15) 



(2.16) 



A) p = is a fixed point of P. 

B) For small \p\, the Poincare map has the following form 

P(p)=Ap + 0(\p\ 2 ), A = X(l), 
where X(t) is the principal matrix solution of the homogeneous system 



p = R(t)p (cf. (2.9) and (2.11)) 



(2.17) 



(2.18) 



C) If the moduli all eigenvalues of A, pi, p,%, ■ ■ ■ , Pn-i, are less than 1, then periodic orbit O is asymptoti- 
cally orbitally stable. 



Next, we turn to the randomly perturbed problem (2.5 ). Let x(t) and x t := {@t,Pt) denote the solutions 



of the deterministic and randomly perturbed problems (2.1 ) and d23j) respectively starting with with initial 



condition xo = (0, p$) £ So,<5i- Using the large deviation estimates for solutions of (2.5 I (cf. Lemma 2.1 of 
Chapter 4 in [18 (3, we have 



s a. ( sup \x t - x(t)\ > S 2 ) < exp{-ca 2 }, 

*6[0,2] 



(2.19) 



for certain constant c > 0, provided a > is sufficiently small. Therefore, with exponentially close to 1 

(2.20) 



probability, the trajectory of (2.5 1 xt returns to S , o,2<5 2 a f ter tne fi rst return time 

T a = iu£{t > 0.5 : 9 t = 1}. 



The first return map p = P(p) is defined by 



P(p) == Pr a 



(2.21) 



The following lemma provides the asymptotic description of the PM. 



Lemma 2.1 in 1181 is stated for a system of the form i 2.5 i with P(x) = id. However, the argument used in the proof of this 



lemma applies to systems with P(x) satisfying the assumptions in i 2. 1 after a suitable modification of the action functional 
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Lemma 2.3. The first return map has the following form 

p = A(I + aB()p + a V + 0(a 2 ,\p\ 2 ), 

where 

A = X(1), B = X- 1 (l)X(l), 

C := £ - &(1) T ?7, £ = / h(s) T dW s , V = I X(t)X- 1 (s)H(s) T dW Sl 
J J 



and 



b(t) T 



a(s) T X(s)X- 1 {t)ds. 



(2.22) 

(2.23) 
(2.24) 

(2.25) 



Remark 2.4. The principal matrix solution of the system with periodic coefficients (2.18 1 can be written 
as X(t) = Q(t) ex.p{tA}, where Q(t) is a 1— periodic matrix. Therefore, B in (2.23 1 can be rewritten as 
follows 

B = Q- 1 (0)Q(0) + Q~ 1 {0)AQ{0). (2.26) 



2.4 Stability of the randomly perturbed PM 



If the moduli of all eigenvalues of A (cf. (2.16l), fii,^2, • • • , Pd, are less than 1, then the periodic orbit O 



of the deterministic system (2.1 ) is asymptotically orbitally stable. This means that if the initial condition 
x(0) is chosen sufficiently close to O, the trajectory {x(t),t > 0} will remain the vicinity of O for all 
future times. Even for very small noise intensity a > 0, for any initial condition a generic trajectory of the 



randomly perturbed system (2.5 ) eventually leaves any neighborhood of O due to the large deviations (cf. 
[18]). Nonetheless on finite interval of time t € [0,T], xt with high probability exhibits stable behavior, 
provided a > is sufficiently small and A is a stable matrix. Thus, we expect that the trajectory of the PM 
remains close to the origin for a long time. To describe stability properties of the randomly perturbed PM, 
for the trajectory of (2.22 1 {p n } starting in a measurable set D C M. d containing the origin, we define the 
first exit time 

T( Pn , D) = min{n : p n <£ D}. (2.27) 



In the remainder of this subsection, we formulate two theorems characterizing the distribution of the 



first exit time of the trajectories of the linearized PM (Theorem 2.6 1 and those of the full nonlinear map 



(Theorem 2.7 1. These theorems use certain auxiliary notation, which we review next. 



Definition 2.5. H25V Let Y be a random variable with values in the set of positive integers and let < p < 1. 
We say that Y is asymptotically geometric with parameter p if 



(Y 



n) 



lim 

n^oo P(Y > n) 



(2.28) 



Recall the expression of the PM in ( 2.22 ». Along with the full nonlinear PM (2.22), we also consider its 
linearization 

Q n = A(I + aB( n )g n ^ 1 + ar] n , (2.29) 
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where and {r] n } are Gaussian random processes denned in (2.24 1. 
Suppose A is a stable matrix with the spectral radius 

p(A) <l-e 

for some < e < 1. Then there exists matrix norm in M. dxd such that 

Mil' < l-e, 

(see, for example, [26]). We will refer to the norm in (|2.31[) as the adapted norm. 



Let | • |' denote a vector norm in M compatible with (2.31 1 

\Ax\' < (1 - e)\x\' Vx £ K d , 

and 72 > 7i > such that 



For fixed h > we define 



Ji\x\ < \x\' < 72|x| \/x 6 



D h = {x £ R d : \x\' < h}. 



(2.30) 



(2.31) 



(2.32) 
(2.33) 

(2.34) 



It is instructive first to understand the statistics of the first exit times for the linearized PM (2.29 1 



Theorem 2.6. Suppose (2.301 holds and {g n } denotes the trajectory of (2.291 starting from initial condition 



Qo = 0. Then for certain = O(e) and for all < a < gq, r(g n , D^) is an asymptotically geometric RV 
with parameter 



Cicrexp 



-Co 



[l + 0{a 2 )) <p< C 3 aexp 



-Ca 



a- 



(i + o(- 2 )) 



(2.35) 



for some positive constants Ci 2 3 4 independent from a. 



We are now in a position to formulate the main result of this paper. 



Theorem 2.7. Suppose (2.30) holds and {p n } denotes the trajectory of (2.22 1 starting from initial condition 



po = 0. Then for certain o~o = 0(e) and for all < a < o~q, RV T(g n , Dh) is subject to the following 
estimate: for n > 1, 

P(f = n + 1) 



P(f >n + 1) 
where C is a positive constant independent of a and h. 



< p, where 0<p< C(ah 2 ) 2/3 , 



(2.36) 



3 The derivation of the Poincare map 

In this section, we compute the linear part of the Poincare map of the periodic orbit O. 
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3.1 The variational equation 



The first step in the derivation of the PM is the asymptotic approximation of the solution of (3.6 1 and (3.7 1 
subject to initial condition 

(9 = and \p \ < 6. (3.1) 

Lemma 3.1. On any finite interval of time, for sufficiently small a > 0, the solution of the initial value 
problem ( |J.6| ), < \3. 7| ), and 1\ admits the following asymptotic expansion 

t = t- b(t) T X(t)p + o-(Ct + O(\po\)) + 0(a 2 , \ P0 \ 2 ), and p t = X(t)p + am + 0(a 2 , \p \ 2 ), (3.2) 

where ^ ^ 

ru= I X(t,s)H(s) J dW s , 6= / h(s) J dW s , and ( t = & - b(t) T Vu (3.3) 
Jo Jo 



and 



b(8) T = - / a T (s)X(s)X- 1 (e)ds. 
Jo 



Proof: For integration of (3.6 1 and (3.7 ), it is convenient to change 9 to a new variable 

(j) = 6 + b(0) T p, 



where b(t) is defined in (3.4 1. In new coordinates (2. 10 1 and ( |2.1 1 ) become 

it = l + 0(\p t \ 2 )+o-(h(cf> t ) r + 0(\p t \)) W u 
p t = R(c/ ) t)pt + 0(\pt\ 2 ) + a^H( ( p t ) T + 0(\p t \))w t . 
The corresponding initial condition is 

4>o = and |po| < <$• 
On a finite time interval t G [0, 2], we expand xt = (<pt, pt) in the asymptotic sum 



(i) 



(3.4) 

(3.5) 

(3.6) 
(3-7) 

(3.8) 

(3.9) 



where z. 



(0,1) 



[<fit^ > ) • Function zf^ is the deterministic. The first order correction z*jp is a Gaussian 
process. Below we state the corresponding initial value problems for z^. ' 1 ^. The remainder satisfies the 
following estimate (cf. Theorem 2.2, |[T8l ) 



E{ sup \1Z{t,a)\ 2 } < Co 2 
*e[o,2] 



(3.10) 



for some positive constant C. We denote the asymptotic relation in (3.10) by TZ(t, a) = 0{a 2 ). 
The zeroth order problem is given by 



^ 0) = l + O(M 2 ) and pf ] =R(<f>)pf ) +0(\p t \ 2 ), 



(3.11) 



subject to = (0, po) T ■ By integrating (3.11 



, we have 

t + O(\p \ 2 ), and pf ] = X(t)po + O(\p 



The first order problem has the following form 

^ = h(t) T W t and $ = RWpP + H(t) T W t , 



with initial condition x t = (0, 0) T . From (3.13 1, we find 



(3.12) 



(3.13) 



Z t :=4 1)= f Hs) T dW s and r, t := pf ] = [ X{t, s)H(s) 1 dW s . (3.14) 
Jo Jo 



Combining ( 3.9 1, ( 3. 12 1 and ( 3. 14 1 and switching back to the original variables (9 t , p t ) (cf. ( 3.5 1), we obtain 

□ 



3.2 Proof of Lemma 1231 



First, we estimate the time of the first return. 
Lemma 3.2. The first return time is given by 

r a = l + 6(l) T X(l) Po - aCi + o(a) + O(\p \ 2 ), 



(3.15) 



Proof: Using the definition of the first return time ( 2.20[ ) and the asymptotic expansion of 9 t ( 3.2 1, we have 

T*-ai(T a )po + Cr a +0(2) = l, a.s., (3.16) 



where ai(t) := b(t) J X(t) and 0(2) := 0(a 2 , \po\ 2 ). It follows from that r a ->■ r a.s., as a ->■ 0, 
where tq is the first return time for the unperturbed deterministic trajectory (er = 0). Thus, we write 



To = to + n(c), 
where ri(cr) = o(l). By plugging ( |3.17 ) in ( |3.15 ), we have 

to + Tx(a) - ai(r )/3 - a' 1 (t')n(o-)po + °(t + o"(Cr CT - Cr ) + 0(2) = 1, 



(3.17) 



(3.18) 



where if is a real number lying between tq and r^. By setting a = in (3.18 1, we obtain 



and, thus, 
Next, 



to - ai(T )po + O(|po| 2 ) = 1, 
r = l + a 1 (l)p + 0(|po| 2 ). 

n(«7)(l - fll(t')po) + ffCl + ^{(Cr CT - Cr ) + (Cr " Cl)} = 0. 



(3.19) 



By the Ito isometry, two terms in the curly brackets are o(l) and O(|po|) respectively. Thus, 

n(a) = -aCi + O{a\p \) + o(a). 



The combination of (3.17 1, (3.19), and ( |3.20 i proves the lemma. 

□ 



(3.20) 



We are now in a position to compute the first return map. From ( |3.2[ ) and ( |3.15| >, we have 
p = X(T a ) Po + ar, Ta + 0(2) = X(l + ctCi)po + <"/i + 0(2). 

Further, 

X(l + ad) = X(1)(I + aX-^XilXi) + 0(a 2 ). 



This proves (2.22). 



(3.21) 



4 The randomly perturbed linear map 



In this section, we study stability of the linearized PM (2.29 1. Theorem 2.6 will follow from the analysis of 
the following slightly more general class of equations 



where 



y n = M n y n -i + z n , n> 1, 



A(I + ai n B) and z n = 5Grj n . 



(4.1) 



(4.2) 



Here, A, B,G E R dxd ; £ n and r] n are independent identically distributed (IID) copies of standard normal 
RV £ and rj in E and ~R d respectively. Furthermore, we assume that A is a stable matrix with spectral radius 



subject to (2.30 1 and G is nondegenerate. In fact, without loss of generality, we assume that G is positive 
definite, since otherwise, one can take G := (GG T ) 1//2 without changing the distribution of z n . 



Theorem 4.1. Let h > be fixed. For yo E D^, denote the first exit time of the trajectory of (4.1 1 from 
Dh (cf. {2.34)) by r := r(y n , D^). Then there exist positive Sq = 0(e) and uq = 0(e) such that for any 



yo E Dh, o~ E (0, (To), arcc? 5 E (0,5o), the first exit time r := T{y n ,Dh) is asymptotically geometric RV 
with parameter p satisfying 



Cicrexp 



-Co 



;i + 0(a 2 )) <p<C 3 aexp 



-Ca 



a- 



(i + o(- 2 )) 



(4.3) 



for some positive constants C% 2 3 4, which do not depend on a and 5. 



Remark 4.2. Note that setting 5 := a and with appropriate choice of G, (4.1 1 coincides with the expresion 



of the linearized PM (2.29 1. Therefore, Theorem 2.6 is a particular case of Theorem 4.1 



The proof of Theorem |4. 1 1 relies on the convergence results for iterative processes due to Kesten [28]. 
Specifically, our assumptions on {M n } and {z n } guarantee that {y n } converge in distribution to a random 
vector y E M. d as n — > 00. Furthermore, £ and z are independent of y. We verify the conditions of the 



Kesten's theorem implying the convergence of {y n } in Lemma 4.3 
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Lemma 4.3. Under our assumptions on the coefficients in (4.1 ), stochastic process {y n } converges in dis- 
tribution to a random vector y G M. d as n — > oo, such that (£, z) are independent of y. 



The gist of the proof of Theorem 4.1 lies in estimating ¥(My + z ^ This is the subject of the 

following lemma. 



Lemma 4.4. There exist constants < C5 < Cg < 1 swc/z f/iaf, uniformly over s G -Dfr 

C 5 < P(Ms + z £ L> A ) < C 6 . 



(4.4) 



We first prove Theorem 4.1 using Lemma 4.4 and then give the more technical proof of Lemma 4.4 



(r = n) 



followed by the proof of auxiliary Lemma 4.3 
Proof: (Theorem |4.1| ) We need to show 

lim 

n^oo P(r > n) 
and estimate p in terms of the coefficients of Equation ( |4~T 
First, we note 

P( T = n + l) = 



p > 



n+i $ D h ,y k e D h ,k e [n]) 
n+1 i D h \y k G D h ,k £ [n])F(y k G D^, k G [n]) 
1 <£D h \y k eD h ,ke [n])P(r>n + l). 

Further, using ( |4.1| ) and ( |4.2[ ) we have 

P(y n+ i £ D fe |y fc €D h ,k£ [n]) = F(M n+1 y n + z n D h \y n G D ft ). 
The combination of ( |4.6| > and ( |4.7| ) shows that ( |4.5[ ) is equivalent to convergence of 

P(M n+1 y n + z n i D h , y n G D h ) 



¥{M n+l y n + z n i D h \y n G D h ) 



n G D h ) 



to a nonzero limit. 



By Lemma 4.3 



(yn,tn,Zn) -> (y,(,,z), n-^ OO, 



where (£, z) is independent of y. Thus, the limit of the numerator in (4.8 ) as n — > 00 is 



Jd h 



where F y (-) stands for the distribution function of y. By Lemma 



4.4 



C 5 F{yeD h )< [ F(Ms + z^D h )dF y (s)<C e F(y€D h ), 



(4.5) 



(4.6) 



(4.7) 



(4.8) 



(4.9) 
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and, therefore, the limit, p 



= F(My + z £ D h \y G D h ), on the right-hand side of ( 148] ) satisfies 

< C 5 < p < C 6 < 1. 



□ 



Proof: (Lemma [44]) 



1. Let s G Z^, M' < fo. Then, using (2.32 1 and (4.2), we have 



PflMs + z\' >h)< F(a\£\\ABs\' + 6\Gtj\' > eh). 
Recall that £ € A/"(0, 1) and define event 

'\Z\< 



J 7 



2ct AB ' 



Using (4.11 1, we take the estimate in (4.10l one step further 



|.V«-r :|'> />) <P( IGr/l' > + 



(4.10) 



(4.11) 



(4.12) 



We bound the second term on the right hand side of (4.12), using the normal distribution of £ G 

A/-(0,1) 



2 a 
< a/-— exp 

7T GfC 



-C 2 e 2 
2cj 2 



;i + 0(-V 1 )), 



(4.13) 



where C7 = (2|| AB||') is independent from a. 



2. Our next goal is to bound the first term on the right hand side of (4. 12 1. Gaussian random vector 

z = Gr] (4.14) 



has zero mean and co variance matrix GG T . Denote the median of the distribution of z by m and 



K = {x G R d : \Gx\ < m} and K r = {x G R d : dist (x, K) < r}. 



(4.15) 



We will need the following concentration inequality for Gaussian random vectors (cf. QUI pp. 19-21]) 

exp{— 2 _1 r 2 } 



[z G K c ;) < 1 - $(r) 



'2irr 



(l + 0(r~ 2 )), 



(4.16) 



where <£(r) = ^= J'^ e s2 / 2 ds. We also employ an observation of Kwapien ||29l to bound the 
median of z 



m 2 < E \Gr]\ 2 = E Tr{(Gr]) T (Gr])} = Tr E (G W T G T ) = 



(4.17) 
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where \\G\\f = y/Tr (GG T ) is the Frobenius norm of G. From (4.17 1 and the triangle inequality, for 
any y £ K r there exists x E K such that 



\Gy\ < \Gx\ + \G(y-x)\ < \\G\\ F + \\G\\r. 



(4.18) 



We are now in a position to use the concentration inequality (4. 16 1 to bound the first term on the right 
hand side of (14. 121). 



First, using (2.33 1 we switch to the Euclidean vector norm 



2 I >26 £ 



eh \ 



Next, we specify the upper bound on S 



0<5< 



eh 



4\\G\\ F j 2 



where 72 is given by ( |2.33| ), and choose 



Then, ( |4.20| ) implies that 



eh 



and (4.21 1 is 



Hence, 



G\ f < 



eh 



U^ 2 \\G\\ 
eh 



(4.19) 



(4.20) 



(4.21) 



4^72 



45 72 ' 



IGllr. 



( h e h ^ ,.„,. 
+ tt — > \\G\\f + \\G\\r. 



2^72 4^72 4(572 

Thus, using (4. 16 1 and (4.22), we obtain that the right-hand side of ( |4.19| ) is bounded by 



(4.22) 



eh \ 



7 E (|:| > ||(7||/.-+ !|(V|k) 7 PUv;.') < 6XP{ J= ^ (l + 0(r- 2 )) . (4.23) 



'27rr 



Finally, by combining (4.19l and (4.23), recalling the definition of z in (4.14), we arrive at 



(4.24) 



where r is specified in (4.21 ). 
3. Combining ( |4~T2| ), ( [403] ), and ( |4~24| ), we have 



P(Ms + ^fl/,) < y-— exp 

7T 676 



2 (T 



+ 



2vrC 8 e/i 



cxp 



-C 7 2 e 2 
2a 2 

2<5 2 



l + O^V 1 )) 

1(1 + 0(^-2^)) 
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4. To get a lower bound on F(Ms + z ^ Dh), we invoke the Anderson inequality HI 

P ((z - a?) € £>ft) < P(* G Vx e M d , (4.25) 
which is applicable since is convex and symmetric about 0. By ( |4.25 ), 



F(Ms + z i D h ) > P(|Cr?| > toT 1 ). (4.26) 

Further, 

P(|Gty| > M" 1 ) > P(|??| > /i(Ai(G)<5)- x ) > 

P^I^MAi^)- 1 ) = f-^v{^-){l + 0{5^)). (4.27) 

□ 



Proof: (Lemma 4.3 1 Recall that according to Kesten's results (see the beginning of Section 3 or Theorem 6 



in [28 ]) for the convergence in distribution of {y n } we need to know that 

n 

E|Gt/| 7 < oo for some 7 > 0, and that a := lim II IT A(I + cr^ k B)\\ < 0. 



n— >oo 

k=l 



(The existence of such a is guaranteed by the condition E log + \\A(I + a£B) \\ < 00, see GDI.) The first of 
these conditions is obviously true and holds for every 7 > 0. To verify the second condition, we next show 



1 n 

lim - log || TT A(I + a£ k B)\\ < 0, (4.28) 

n— >oo n J- 

provided 



a := 

n— >oo 7), 

fc=l 



a < t~~ r • (4.29) 
c(l-e) 



for a constant c > specified below. 



As we mentioned, the existence of a is guaranteed by the integrability of log ||.A(I+£.B) ||. Furthermore, 
by GUI Theorems 1 and 2] it suffices to show that 

1 - 

lim -E log || TT A(I + a£ k B) \\ < 0. (4.30) 

fe=l 

Obviously, we can replace || • || by any other matrix norm since neither condition is affected by passing to 
an equivalent norm. 



We will show that (430 1 holds for the norm || • ||'. Since x — > logx is concave, by Jensen's inequality 



Elog || H A(I + o-£ k B)\\' < logE|| H A(I + a^ k B)\\'. 



k=l k=l 
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Using corresponding of the adapted norm, independence of ^'s and triangle inequality we see that 

n n 

E|| l[A(I + <rt k B)\\' < HE\\A(I + oZ k B)\\' < (\\A\\'(1 + ac)) n , 



k=l k=l 



where c= \\B\\'E\£\ = \\B\\'^/2/^. 

If a > satisfies ( |4.29| > then ( |2.31| ) implies that ||^H'(1 + ac) < 1 for some c> 0. Thus, by ( |4.30| ) the 
(|4T28T> holds with a < log((l - ei)((l + ere)) < 0. □ 



5 The nonlinear estimates 



In this section, we prove Theorem 2.7 To this end, we recall certain notation used in the previous sec- 
tions. The full and linearized PMs (cf. (2.22 1 and ( 2.29[ >, respectively) are given by the following difference 
equations: 



Pn = 
Qn = 

where r(p) = 0(\p\ 2 ,a 2 ). 
For fixed h > 0, we define 

and the fist exit time f := r(p n , Dh). 



A(I + aB( n )p n -i + ar] n + r(p n -x) 
A(I + aB(n) 



D h = {x G R d : \x\' < h}. 



Our goal is to prove the statement of the Theorem 2.7 for n > 1, 

P(r = n+ 1) 



(f > n+ 1) 



<p, where < p < C(ah 2 ) 2/3 . 



(5.1) 
(5.2) 



(5.3) 



As in the proof of Theorem 4. 1 we obtain 



(5.4) 



Denote x n = p n — g n and consider the difference equation 

A(I + a^ n B)x n -i + r(p n _i] 



Iterating (5.5 1, we obtain 



/ n— 1 



n+1 /n-i 



Vfc=0 / j=2 Vfc=0 / 



(5.5) 



(5.6) 
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where we used the convention on the right-hand side that the product over the empty range is 1. Note that 
n < f implies that for j = 2, . . . , n + 1 we have \pj-2\' < h. Hence, there exists a constant C such that for 
all such j's |r(yj_2)|' < Ch 2 . Also, since xq = the first term in ( |5.6| > vanishes. Furthermore, for t € M. 

and w G M. d 

\A(I + tB)w\' < \Aw\' + \A(tBw)\' < (1 - e)(|^|' + \tBw\') < (1 - e)(l + \t\ • \\B\\')\w\\ 
which implies that 



n-j 



Y[ A(I + aU-kB)r(p 



fc=0 



n-j 



< a - na + • \\b\\'m Pj _ 2 )\>. 



k=0 



Therefore, whenever n < r, we have 



r < ch 2 £(i - e )^' +i n(i + 'ft,.*! ■ w) < ^ 2 E( x - e i)' + • W)- 

j=2 k=0 1=0 k=l 



where C = C(e, cr, A, B) is an absolute constant whose value may vary from use to use. We bound the 
variance of the latter sum as follows: 



var £(1 - eYl[(l + a\£ k \ ■ \\B\\ f ) = - e)^ var JJ(1 + • 

\£=0 fc=l / £=1 \fc=l 

oo oo / ^ m \ 

+ 2 E E (i-^ +m co V na+^i-iiBin.n^+^i-Pin 

£=1 m=£+l \k=l k=l J 

oo oo / £ m \ 

^ 2 E X> - ^ +m cov na + ■ \\ B n + ■ iw) • 



fcl m =£ \ 

By independence, for m > I > we have 

£ m 



k=l 



k=l 



e n(! + • ii fl ir) + • 11*10 = w + • iisiiwu + ^r-', (5.7) 



*;=i 

where we have set 



k=l 



fi B = \\B\\'-E\t\ = \\B\\'^2/tt. 
Therefore, the covariance above is equal to 

(E(l + a\Z\ ■ \\B\\') 2 Y(1 + ansr-' ~ (1 + <?VbT +1 

= (1 + o^bY 1 - 1 {(E(l + a|f | • USH') 2 / - ((1 + W)'} . 

The term in the curly brackets is bounded above by 

(E(l + ■ \\B\\) 2 - (1 + a„ B ) 2 ) £(E(1 + <r|£| ■ ||S||') 2 y((l + ^ B f) l ~^ . 

j=0 



(5.8) 



(5.9) 
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The factor outside the summation is equal to var(cr||i?|| / • |£|) = 0(a 2 ). To bound the sum in (5^9) we use 
the inequalities 

(1 + ap B ) 2 < E(l + o-|e| • ||5||0 2 = 1 + 2a\\B\\'^ + {a\\B\\'f < (1 + a||B||') 2 
which imply that the sum above is bounded by 

e-i 

2(1 + °\\B\\') 2j (l + aWBW') 2 ^- 1 -^ < i{\ + aWBW') 2 ^. 

j=0 



Combining ( |5.7| ) and ( |5.9| ) with the expression bounding the variance we get 

oo oo 

Co 2 " e) 2l (l + a\\B\\r^ £((1 ~ e )U + ^)) 



m—i 



(5.10) 



m=i 



Using ( |5.8| > we see that if er < oo where, say, do < 6/211511' then the double sum in ( 5- 10 > is bounded by a 
constant depending on e only, and thus 

/ oo £ \ 

var £(l-e)'n(l + a|&|. WO < ^ 

V=0 fc=l / 

Therefore, by Chebyshev's inequality, for any hi > we have 

n+l n-j 



3'/ 



J=2 



fc=0 



(5.11) 



Note that if hi = 0{a a h^) for < a < 1 and 1 < /3 < 3/2. (The optimal values of a and j3 will be 
identified below.) Then the right hand side of (5.11 ) is 0(a 2<yl ~ a ^ h 2 ( 2 ~^). For such hi, set h = h — h\. 
Since p n = g n + x n we have 

\'<h) = 



+ ar„|' <h)> P(|y n |' <h—h%, \x n \' < hi) 
Vn\' < h )-F(\y n \' < ho,\xn\' > hi) 
> P(|y«r < ho) - P(\x n \' > hi) 



This bounds from below the denominator in ( |5.4[ ). To upper bound the numerator P(|p n+ i |' > h,\p n \' < h) 
we first write 



'n+l 



' > h, \ Pn \' <h)< F{\y n+ i\' >h-hi,\y n + x n \' < h) + P(|x n+1 |' > hi). (5.12) 



Now, using that \y n \' — \x n \' < \y n + x n \' we see that \y n + x n \ < h implies that either \x n \' > hi or 
\yn\' < h + hi. Therefore, 

P(|y n+ i|' >h-hi,\y n + x n \' < h)<F(\y n+ i\' > h ,\y n \' < h + hi) + F(\x n \' > hi) 
< F(\y n+ i\' > h , \y n \' < h ) + F(h - hi < \y n \' <h + hi) + F(\x n \' > hi) 
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Combining this with ( |5.12[ ) we obtain 



PdPn+il' > h,\Vn\' < h l K gCk+lT > h ,\y n \' < h )+F{h-h 1 < \y n \' < h + h x ) + 2P(|x n |' > h x ) 



\Pn\' < h) 



\Vn\' < ho)-F(\x n \> > h ± ) 



Using the inequality < ^f 1 valid for all t > 0, < s < t and < u < t we can bound the last 
expression by 

F(\y n+1 \' > hp, \y n \' < h ) + F(h - hi < \y n \' <h + h 1 ) + 3P(|z n |' > hj) 

n\y n \' < ho) 

Furthermore, since y n is a continuous random variable, x —> F(\y n \' < x) is uniformly continuous and since 
hx = 0{o a hP) andP(|y n |' < ho) = 0(1) andP(|x n |' > h) = 0(a 2 ^-^h 2( - 2 ~^) we can bound the last 
quantity by 

F(\y n+1 \'>h ,\y n \'<h ) + Q{aahP) + 0(a2(1 _ Q) ^ M) 



\Vn\' < ho) 



(5.13) 



To optimize the last two terms on the right hand side of (5.13 1, we ch oose a = 2/3 and /3 = 4/3, which 

is of much 



makes them 0((ah 2 ) 2 ^ 3 ). The first term on the right-hand side od (5.13 1, by Theorem 



2.6 



smaller order. This concludes the proof of Theorem 2.7 



6 Applications 



In this section, we illustrate the analysis in the previous sections with applications to three problems in math- 
ematical neuroscience. To this end, we use two conductance-based models of single neurons and a model of 
electrically coupled network of pancreatic /3— cells. In the absence of noise, all these models exhibit stable 
periodic oscillations. Adding noise to these models results in more complex stochastic oscillatory regimes: 
irregular bursting and mixed-mode oscillations. We show that despite different mathematical formulations 
of these models and different forms of resultant oscillations, emergent stochastic oscillatory patterns in these 
models are formed due to the random perturbations of stable limit cycle oscillations. We show numerically 



that the number of spikes in one burst in models in <:6.1 and {6.2 and the number of small oscillations in 



<: 6.3 are distributed approximately geometrically in accord with Theorem 2.7 



Numerical examples in §6.1| and ^2] appeared before in E51I37I respectively. We use them here to give 
the reader a feeling for the range of possible applications of our analytical results. The numerical example 



in £ 6.3 is new. In somewhat different dynamical regime for a stochastically forced 2D FitzHugh-Nagumo 
oscillator, asymptotically geometric distribution of small amplitude oscillations was shown in [9l (see also 
[38 ] for related results). Finally, we note that there is another mechanism for generating irregular bursting 
and mixed-mode oscillations featuring geometric distribution. It does not involve random perturbations, but 
is based on certain properties of chaotic attractors in closely related neuronal models (see P4"ll36l0 . 



6.1 Bursting in a single cell model 



Our first example concerns bursting. This is a common pattern of electrical activity in neurons. Bursting 
is characterized by series of fast oscillations separated by periods of no activity (Fig. [1J>). Understanding 
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(0 500 900 

Number of spikes 



Figure 1: Randomly perturbed model of neural oscillator (6.1 1-(6.3 1, (a) The projection of a periodic tra- 
jectory of the unperturbed deterministic system (red) and that of the randomly perturbed one (blue) onto 
n-v plane, (b) The timeseries generated by the randomly perturbed model: series of fast oscillations closely 
resembling oscillations of the unperturbed deterministic system are separated by periods of no activity, (c) 
The tail of the normalized histogram of the number of spikes in one burst is fitted with an exponential curve. 



dynamical mechanisms underlying bursting is an important problem in mathematical neuroscience ETTl . 
Below, we describe a mechanism of stochastic bursting based on random perturbations of a stable limit 
cycle. 

The following system of three differential equations describes the dynamics of a hypothetical neuron 
(cf.§9.1.1, ED): 

C m v = -gNapmoo{v){v - E NaP ) - g K n(v - E K ) - g K Mv(v - E K ) (6.1) 

- g L (v - E L ) + I + aw t , 
"oo(v)-n 

n = , (6.2) 



Vooiv) - y 



(6.3) 



Here, v, n, and y stand for the cell membrane potential and two gating variables respectively. Param- 
eters g s and E s are the maximal conductance and the reversal potential of different ionic currents I s , 
s G {NaP, K, KM, L} and I is the applied current. The time constants r n and r y determine the rates 
of activation in the populations of K and KM channels. The steady-state functions are defined by s^v) = 

^1 + exp {^Y^j ^ , s G {m, n, y} . In addition, a white noise process of intensity a is added to the right 
hand side of the voltage equation to model various fluctuations affecting membrane potential. 

The parameter values are summarized in the following table. 



Table 



gNa 


20 mS/ cm 2 


gK 


10 mS/cm 2 


gxM 


5 mS/cm 2 


gL 


8 mS/cm 2 


ENa 


60 mV 


Er 


-90 mV 


Ei 


-80 mV 




0.152 ms- 1 




20 ms^ 1 


a 


1 


I 


5pA 




-20 mV 


a n 


-25 mV 


a y 


-10 mV 
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5 


by 


5 


Cm 


1 nF/cm 2 
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Figure 2: Randomly perturbed model of two coupled neural oscillators ( |6-4| )-( [6T6[ >. (a) The projection of a 
periodic trajectory of the unperturbed deterministic system (red) and that of a randomly perturbed one (blue) 
onto n-v plane, (b) The timeseries generated by the randomly perturbed coupled model, (c) The normalized 
histogram of the number of spikes in one synchronized burst features asymptotically geometric distribution. 



The deterministic model (6.1 1-(6.3 1 with a = has a stable limit cycle L (see Fig. [Tp). Under small 



random perturbation (a > 0), the trajectory of (6.1 1-(6.3) after a random number of rotations around the 



limit cycle of the deterministic system leaves the basin of attraction of the limit cycle and undergoes a large 
excursion in the phase space before coming back to a small neighborhood of L. By applying the results of 
this paper to the Poincare map of the limit cycle for this problem, one can show that the number of spikes in 
one burst is distributed approximately geometrically. The tail of the distribution of the number of spikes in 
one burst shown in Fig. [I]: is well fitted by an exponential function. 



Stochastic bursting in conductance-based models similar to (|6. l|>-([673]> was analyzed in |[25l . The results 



of the paper provide a more general treatment of noise-induced bursting. Furthermore, the examples in the 
remainder of this section lie outside the scope of applicability of the method in [25 ], but can be effectively 
analyzed using the results of this paper. 



6.2 Bursting in a coupled network 

Our next example presents a model of a coupled network generating synchronous bursting oscillations. 
This is a model of an ensemble of pancreatic (3— cells coupled electrically. For simplicity, we use a two- 
cell network and refer the interested reader to ||3"7l for the description and numerical simulations of bigger 
networks. 

The dynamics of Cell i is governed by the following system of differential equations 

C m 0«> = -/ io „(»"»,n«>, !/ l i l) + 9 ( t ,«-^«>) + ™l'>, (6.4) 

t(vW) 

y« = e(Ica(v {i) ) - fcy«), i € {1, 2}, j = 3 - i. (6.6) 



The biophysical meaning of these equations is similar to that of the model in £6.1 i>M, rW, and j/W stand 
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for the membrane potential, gating variable, and calcium concentration corresponding to Cell i. Random 
processes i € {1, 2}, are independent copies of standard Brownian motion. 



The first term on the right hand side of ( 6.4 1 models the combined effect of sodium and calcium currents, 
iNa+Ca, the calcium-dependent potassium current, Ixca, delayed rectifier 1% and a small leak current, // 



iNa+Ca + T-KCa + Ir + h, 



(6.7) 



where 



INa+Ca 

Irc 
Ik 



gNa+Cam 00 (v) 3 h 00 (v)(E I - v), 

gKCaU , v 

g K n A (E K - v), 
gi(E t - v). 



The steady state functions used to model the ionic currents above are given by 

fooiv) = T^rl < V f£{ m ,h,n}, 
otf(v) + Pf(v) 



where 

a m = 0.1(u + 25)(l -exp{-0.1(u + 25)})- 1 ,/3 m = 4exp{-(w + 50)/18}, a h = 0.07exp{-0.05(t; + 50)}, 
P h = (1 + exp{-0.1(v + 20)}) _1 ,a n = 0.01(w + 20)(1 - exp{-0.1(u + 20)})" 1 , f3 n = 0.125 exp{-(v + 30)/80}. 

The time constant of the delayed rectifier is given by 

T n = (230(a n + /3 n ))" 1 . 

The values of the remaining parameters are summarized in the following table. 



Table 



9Na+Ca 


1800s- 1 


ENa+Ca 


\WmV 


9K 


1700s" 1 


Er 


-15mV 


kc 


1 [ 8 1 mv 


9l 


Is- 1 


Ei 


-40mV 


9RC 


12s" 1 


Eca 


\00mV 


e 




Cm 


lfiF/cm 2 


a 


10 mVs- 1 



When uncoupled (g = 0), the models of each cell generate stable oscillations similar to those generated 
by the model in j ^6.1| One can show that for sufficiently strong coupling (g 3> 1), the deterministic coupled 
system (a = 0) has a stable limit cycle corresponding to synchronous oscillations in both cells (cf. II3510 . 
In the presence of noise (a > 0), the trajectory of the coupled system once in a while leaves the basin of 
attraction of the limit cycle. This terminates one burst. Since one spike in a burst corespond to one iteration 
of the PM, we conclude that the number of spikes in one synchronized burst of the coupled system must be 
distributed approximatelly geometrically as shown in Fig. [2] 
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6.3 Mixed-mode oscillations 



Our final example deals with mixed-mode oscillations, another type of nonlinear oscillations that are im- 
portant in neuroscience |[T5l l9l [36ll38l . To this end, we use a modification of the Hodgkin-Huxley model of 
a neuron in the regime close to the Andronov-Hopf bifurcation. It was introduced by Doi and Kumagai in 
fT5l . The model consists of the differential equations for the membrane potential v and two gating variables 
n and h: 



C m v 
h 

h 



-gN a moo(v) 3 h(v - E Na ) - g k n 4 (v - E K ) - g t {v - E{) + aw, 

T n (v) 

hoojv) - h 

Th{v) 



(6.8) 
(6.9) 

(6.10) 



The biophysical meaning of these equations is similar to that of the model discussed in { 6. 1 The steady 
state functions 

OLf{v) 



foo(v) = — , s , r , fe{m,h,n}, 
a f {v)+l3f{v) 



and the time constants 



are defined using 



Tf(v) 



a f (y) + (3f(v) 



a m = 0.1(« - 25)(1 - exp{0.1(25 - v)})' 1 , f3 m = 4exp{-t;/18}, a h = 0.07exp{-0.05u}, 
Ph = (1 + exp{0.1(30 -v)})~ l ,a n = 0.01(« - 10)(1 - exp{0.1(10 - v)})^ 1 , f3 n = 0.125 exp{-v/80}. 

The values of the remaining parameters are summarized in the following table. 



Table 



9Na 


120 mS/cm 2 


9K 


36 mS/cm 2 


9l 


0.3mS/cm 2 


ENa 


115 mV 


Er 


-12 mV 


Ei 


10.5999my 


f n 


20 


fh 


1 


a 


5 • 10" 4 







The parameters in the neuronal model (6.8 1-(6.10) are chosen such that the deterministic model (a = 0) 
has a stable limit cycle, whose projection onto v — n plane is shown in Fig. [3^. Under the action of noise, 
after a random number of rotations around the periodic orbit of the deterministic system, the trajectory leaves 
the vicinity of the limit cycle. The global structure of the vector field of the deterministic model guarantees 
that the trajectory of the randomly perturbed system returns to the neighborhood of the limit cycle after each 
excursion in the phase space. This results in mixed-mode oscillations of the membrane potential consisting 
of a random number of small oscillations separated by large spikes (see Fig. [3j)). In accord with the results 
of this paper, we find that the distribution of the number of small oscillations generated by this model is 
approximately geometric (see Fig. [3b). 
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Figure 3: Noise-induced mixed-mode oscillations generated by (6.8 l-(6.10l. (a) The projection of a peri- 
odic trajectory of the unperturbed deterministic system (red) and that of a randomly perturbed one (blue) 
onto n-v plane, (b) The timeseries generated by the randomly perturbed coupled model, (c) The normal- 
ized histogram of the number of small oscillations between two consecutive spikes features asymptotically 
geometric distribution. 



7 Discussion 



Since the time of Lyapunov and Poincare, problems in nonlinear oscillations stimulated and guided the de- 
velopment of the geometric theory of ordinary differential equations and its diverse applications to physics 
and biology |f2j|22l|32l. Many effective analytical methods have been developed for studying periodic mo- 
tion and effects of forcing in deterministic systems (6) El [23]]. In contrast, apart from the large deviation 
type techniques ifTHl and those for slow-fast systems [8], there are few general analytical approaches avail- 
able for studying effects of random forcing on nonlinear oscillations and the analysis of such systems is 
often done on the case-by-case basis |[5ll3ll9l [T4ll25l . 

In the geometric theory of differential equations, the principal tool for studying stability of periodic 
motion is the reduction to a PM. This work presents a systematic construction of the PM for an important 
class of randomly perturbed problems: a limit cycle oscillator forced by small white noise. For trajectories 
of so-obtained PM, we analyzed the statistics of the first exit times from a small neighborhood of the origin, 
which corresponds to the periodic solution of the unperturbed system. We showed that if the periodic 
solution of the deterministic system is asymptotically stable the first exit times have approximately geometric 
distribution. This result implies universality of the geometric distribution in diverse oscillatory regimes 
generated by random perturbations of stable limit cycles. 

Irregular oscillations featuring this dynamical mechanism are common for differential equation mod- 
els in applied science and, in particular, in mathematical biology. We provided three representative ex- 
amples from biophysics: irregular bursting and mixed-mode oscillations generated by randomly perturbed 
conductance-based models of neurons and synchronous noise-induced bursting in randomly forced neuronal 
network. We showed that these dynamical regimes are caused by random perturbations of stable limit cycles 
and, therefore, all of them feature geometric distribution. 
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